clear all;
clear global;
close all;
format short;
echo on;

%% Add the paths for the toolbox
addpath('C:\MATLAB6p5\spatial_econometrics\regress')
addpath('C:\MATLAB6p5\spatial_econometrics\distrib')
addpath('C:\MATLAB6p5\spatial_econometrics\util')
addpath('C:\MATLAB6p5\spatial_econometrics\graphs')
addpath('C:\MATLAB6p5\spatial_econometrics\optimize')
addpath('C:\MATLAB6p5\spatial_econometrics\spatial\sar_models')
addpath('C:\MATLAB6p5\spatial_econometrics\spatial\weights')
addpath('C:\MATLAB6p5\spatial_econometrics\spatial\lndet_funcs')
addpath('C:\MATLAB6p5\spatial_econometrics\spatial\sac_models')
addpath('C:\MATLAB6p5\spatial_econometrics\spatial\saw')
addpath('C:\MATLAB6p5\spatial_econometrics\spatial\gmm')
addpath('C:\MATLAB6p5\spatial_econometrics\spatial\sem_models')
addpath('C:\MATLAB6p5\spatial_econometrics\spatial\stats')
addpath('C:\MATLAB6p5\spatial_econometrics\spatial\newsaw')

%% Add the paths for the work directotry
addpath('C:\MATLAB6p5\work\Beck_etal_archive\')


% Table 1 analysis

% load data and matrix
dat98 = CSVREAD('C:\MATLAB6p5\work\Beck_etal_archive\basedata98.csv',1,0);
dmat98 = CSVREAD('C:\MATLAB6p5\work\Beck_etal_archive\dmat98.csv',1,0);
tmat98 = CSVREAD('C:\MATLAB6p5\work\Beck_etal_archive\tmat98.csv',1,0);

% Make distance a standardized sparse matrix 
bdmat = sparse(dmat98(:,1),dmat98(:,2),dmat98(:,3)); 
rowsum = bdmat * ones(length(bdmat(:,1)),1);
ndmat = spdiags(1./rowsum,0,length(bdmat(:,1)),length(bdmat(:,1)))*bdmat;
ndmat2=ndmat/2;

% Make trade a standardized sparse matrix 
bt3 = sparse(tmat98(:,1),tmat98(:,2),tmat98(:,3)); 
rowsum = bt3 * ones(length(bt3(:,1)),1);
nt3 = spdiags(1./rowsum,0,length(bt3(:,1)),length(bt3(:,1)))*bt3;
nt3_2=nt3/2;

%  create data matrix
y98 = dat98(:,2);
x=dat98(:,3);
xc98 = [ ones(length(dat98(:,1)),1),x];

cvnames = strvcat('dma','constant','loggdppc');
 
%diagnostics
morandistance = moran(y98,xc98,ndmat);
prt(morandistance);
morantrade = moran(y98,xc98,nt3);
prt(morantrade);
lmsardistance=lmsar(y98,xc98,ndmat,ndmat);
prt(lmsardistance);
lmsartrade=lmsar(y98,xc98,nt3,nt3);
prt(lmsartrade);
lmerrordistance=lmerror(y98,xc98,ndmat);
prt(lmerrordistance);
lmerrortrade=lmerror(y98,xc98,nt3);
prt(lmerrortrade);

% ols
ols98 = ols(y98,xc98);
prt(ols98,cvnames);

% sar with distance
distsar98= sar(y98,xc98,ndmat);
prt(distsar98,cvnames);

% distance sar lr effects
a=(inv(speye(170)-(.474989*ndmat)));
changey=zeros(170,1);
for i=1:170,
   changex=[zeros((i-1),1);1;zeros((170-i),1)];
   temp=a*changex*1.725286;
   changey(i)=temp(i);
end
mean(changey) %1.8048
min(changey) %1.7500
max(changey) %2.0373

% sar with trade
tradesar98=sar(y98,xc98,nt3);
prt(tradesar98,cvnames);

% trade sar lr effects
a=(inv(speye(170)-(.503997*nt3)));
changey=zeros(170,1);
for i=1:170,
   changex=[zeros((i-1),1);1;zeros((170-i),1)];
   temp=a*changex*2.224202;
   changey(i)=temp(i);
end
mean(changey) %2.2385
min(changey) %2.2243
max(changey) %2.5233

% saw with trade and distance
tradedistsaw98_2=saw2(y98,xc98,nt3_2,ndmat2);
prt_saw(tradedistsaw98_2,cvnames);

%saw lr effects
a=(inv(speye(170)-(.586811*nt3_2)-(.891655*ndmat2)));
changey=zeros(170,1);
for i=1:170,
   changex=[zeros((i-1),1);1;zeros((170-i),1)];
   temp=a*changex*1.528032;
   changey(i)=temp(i);
end
mean(changey) %1.6262
min(changey) %1.5570
max(changey) %2.0193


%Table 2 analysis

%Load data
cmat98_BA = CSVREAD('C:\MATLAB6p5\work\Beck_etal_archive\dyadmat_europe98_BA-no origins 3+.csv',1,0);
cmat98_noBA = CSVREAD('C:\MATLAB6p5\work\Beck_etal_archive\dyadmat_europe98_no BA-no origins 3+.csv',1,0);
logdat98 = CSVREAD('C:\MATLAB6p5\work\Beck_etal_archive\dyaddata_europe98_log-no origins 3+.csv',1,0);
logdat98 = sortrows(logdat98,1);

%Make weight matrices
bcmat_BA = sparse(cmat98_BA(:,1),cmat98_BA(:,2),cmat98_BA(:,3)); 
rowsum = ones(length(bcmat_BA(:,1)),1);
rowsum = bcmat_BA * ones(length(bcmat_BA(:,1)),1);
ncmat_BA = spdiags(1./rowsum,0,length(bcmat_BA(:,1)),length(bcmat_BA(:,1)))*bcmat_BA;
ncmat_BA_2=ncmat_BA/2;

bcmat_noBA = sparse(cmat98_noBA(:,1),cmat98_noBA(:,2),cmat98_noBA(:,3)); 
rowsum = ones(length(bcmat_noBA(:,1)),1);
rowsum = bcmat_noBA * ones(length(bcmat_noBA(:,1)),1);
ncmat_noBA = spdiags(1./rowsum,0,length(bcmat_noBA(:,1)),length(bcmat_noBA(:,1)))*bcmat_noBA;
ncmat_noBA_2=ncmat_noBA/2;

%Make data
y98 = logdat98(:,4);
xc98 = [ ones(length(logdat98(:,1)),1),logdat98(:,5),logdat98(:,6),logdat98(:,7),logdat98(:,8),logdat98(:,9),logdat98(:,10),logdat98(:,11),logdat98(:,12)];

cvnames = strvcat('logtrade','constant','logdem','logapop','logbpop','logargdppc','logbrgdppc','logs','logdist','logmid');

%ols
ols_europe98_gravity = ols(y98,xc98);
prt(ols_europe98_gravity,cvnames);

% sem with two matrices, rolled into W, with .32 as the optimal value.
%diary('C:\MATLAB6p5\work\Beck_etal_archive\sem optimize-europe.log');
%for I=([20:40]/100),
%    W=(I*ncmat_BA)+((1-I)*ncmat_noBA);
%    tempsem=sem(y98,xc98,W);
%    prt(tempsem,cvnames);
%end
%diary close
%%I=.32 is optimal

a=.32;
W=(a*ncmat_BA)+((1-a)*ncmat_noBA);
newsem=sem(y98,xc98,W);
prt(newsem,cvnames);


%Table 3 analysis

%Load data
cmat98_BA = CSVREAD('C:\MATLAB6p5\work\Beck_etal_archive\dyadmat_africa98_BA.csv',1,0);
cmat98_noBA = CSVREAD('C:\MATLAB6p5\work\Beck_etal_archive\dyadmat_africa98_no BA.csv',1,0);
logdat98 = CSVREAD('C:\MATLAB6p5\work\Beck_etal_archive\dyaddata_africa98_log.csv',1,0);
logdat98 = sortrows(logdat98,1);

%Make weight matrices
bcmat_BA = sparse(cmat98_BA(:,1),cmat98_BA(:,2),cmat98_BA(:,3)); 
rowsum1 = bcmat_BA * ones(length(bcmat_BA(:,1)),1);
ncmat_BA = spdiags(1./rowsum1,0,length(bcmat_BA(:,1)),length(bcmat_BA(:,1)))*bcmat_BA;

bcmat_noBA = sparse(cmat98_noBA(:,1),cmat98_noBA(:,2),cmat98_noBA(:,3)); 
rowsum2 = bcmat_noBA * ones(length(bcmat_noBA(:,1)),1);
ncmat_noBA = spdiags(1./rowsum2,0,length(bcmat_noBA(:,1)),length(bcmat_noBA(:,1)))*bcmat_noBA;

%Make data
y98 = logdat98(:,4);
xc98 = [ ones(length(logdat98(:,1)),1),logdat98(:,5),logdat98(:,6),logdat98(:,7),logdat98(:,8),logdat98(:,9),logdat98(:,10),logdat98(:,11),logdat98(:,12)];

cvnames = strvcat('logtrade_1','constant','logdem','logapop','logbpop','logargdppc','logbrgdppc','logs','logdist','logmid');

%ols
ols_africa98_gravity = ols(y98,xc98);
prt(ols_africa98_gravity,cvnames);

%sem with two weights

%find optimal a
%diary('C:\MATLAB6p5\work\Beck_etal_archive\sem optimize.log');
%for I=([30:50]/100),
%    W=(I*ncmat_BA)+((1-I)*ncmat_noBA);
%    tempsem=sem(y98,xc98,W);
%    prt(tempsem,cvnames);
%end
%diary close
%%I=0.42 is optimal

a=.42;
W=(a*ncmat_BA)+((1-a)*ncmat_noBA);
newsem=sem(y98,xc98,W);
prt(newsem,cvnames);


%Table 4 analysis

% load data
dat = CSVREAD('C:\MATLAB6p5\work\Beck_etal_archive\dyaddata_morrow.csv',1,0);
dat = sortrows(dat,1);
cmat = CSVREAD('C:\MATLAB6p5\work\Beck_etal_archive\dyadmat_morrow_1.csv',1,0);

% Form weight matrix
bcmat = sparse(cmat(:,1),cmat(:,2),cmat(:,3)); 
rowsum = ones(length(bcmat(:,1)),1);
rowsum = bcmat * ones(length(bcmat(:,1)),1);
ncmat = spdiags(1./rowsum,0,length(bcmat(:,1)),length(bcmat(:,1)))*bcmat;

%  create data matrix
y = dat(:,5);
lagy=dat(:,16);
wlagy=ncmat*lagy
wy=ncmat*y

%models
xc = [ ones(length(dat(:,1)),1),dat(:,6),dat(:,7),dat(:,8),dat(:,9),dat(:,10),dat(:,11),dat(:,12),dat(:,13),dat(:,14),dat(:,15),dat(:,16)];
cvnames=strvcat('xijl','constant','gnpil','gnpjl','popil','popjl','distancel','taul','demdl','midl','multiall','biall','lagxijl');
ols_morrow_lagdv=ols(y,xc);
prt(ols_morrow_lagdv,cvnames);

xc = [ ones(length(dat(:,1)),1),dat(:,6),dat(:,7),dat(:,8),dat(:,9),dat(:,10),dat(:,11),dat(:,12),dat(:,13),dat(:,14),dat(:,15),dat(:,16),wlagy];
cvnames=strvcat('xijl','constant','gnpil','gnpjl','popil','popjl','distancel','taul','demdl','midl','multiall','biall','lagxijl','Wlagxijl');
ols_morrow_lagdv_wlagdv=ols(y,xc);
prt(ols_morrow_lagdv_wlagdv,cvnames);